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ABSTRACT 

Context. Predicting the emerging X-ray spectra in several astrophysical objects is of great importance, in particular when the 
^ observational data are compared with theoretical models. This requires developing numerical routines for the solution of the 
I ' radiative transfer equation according to the expected physical conditions of the systems under study. 
Aims. We have developed an algorithm solving the radiative transfer equation in the Fokker-Planck approximation when both 
thermal and bulk Comptonization take place. The algorithm is essentially a relaxation method, where stable solutions are 
obtained when the system has reached its steady-state equilibrium. 

Methods. We obtained the solution of the radiative transfer equation in the two-dimensional domain defined by the photon 
energy E and optical depth of the system r using finite-differences for the partial derivatives, and imposing specific boundary 
conditions for the solutions. We treated the case of cylindrical accretion onto a magnetized neutron star. 

Results. We considered a blackbody seed spectrum of photons with exponential distribution across the accretion column and for 
an accretion where the velocity reaches its maximum at the stellar surface and at the top of the accretion column, respectively. 
In both cases higher values of the electron temperature and of the optical depth r produce flatter and harder spectra. Other 
parameters contributing to the spectral formation are the steepness of the vertical velocity profile, the albedo at the star 
surface, and the radius of the accretion column. The latter parameter modifies the emerging spectra in a specular way for the 
two assumed accretion profiles. 

? Conclusions. The algorithm has been implemented in the xspec package for X-ray spectral fitting and is specifically dedicated 
r> ' to the physical framework of accretion at the polar cap of a neutron star with a high magnetic field 10^^ G). This latter case 
^ . is expected to be typical of accreting systems such as X-ray pulsars and supergiant fast X-ray transients. 

Key words. Methods: numerical - X-rays: binaries - Radiative transfer - Magnetic fields 
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1. Introduction 

The solution of the radiative transfer equation (RTE) 
that describes the modification of a seed photon spectrum 
due to Comptonization in a plasma is a much debated 
mathematical problem. The equation in its full form is 
indeed integro-differential (Pomraning 1973) and allows 
for analytical solutions under some particular assump- 
tions, such as electron temperature To = (Titarchuk & 
Zannias 1998) or in the energy domain when the emerg- 
ing spectrum is a powerlaw (Titarchuk & Lyubarskij 
1995). If the photon energy exchange for scattering is low 
(Av/v <IC 1), it is possible to perform a Taylor expan- 
sion of the Comptonization operator around the photon 



initial energy, which transforms the RTE from integro- 
differential to purely differential (Rybicki & Lightman 
1979); this is known as the Fokker-Planck (FP) approx- 
imation. The necessary conditions to allow this mathe- 
matical approach are that the Compton-scattering process 
occurs below the Klein-Nishina regime, namely when the 
electron temperature kT^ is subrelativistic 100 keV), 
and that the optical depth of the Comptonization region 
is T ^ 1. The regime of low temperature and high op- 
tical depth of the plasma indeed ensures that the spec- 
trum is almost isotropized, so that it is possible to use 
the Eddington approximation for the specific intensity 
= J{v) -\- 3V- F(i/), where J and F are the zero and 
first moment of the intensity field, respectively. 
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Moreover, if the plasma is not static but subject to 
dynamical (bulk) motion with velocity v(t), then another 
condition for using the FP approximation is that v{t) 
must be subrelativistic. When all the above restrictions 
are considered, one obtains by computing the first two 
moments of the RTE a main equation that describes the 
shape of the angle-averaged emerging intensity J{v) of 
the Comptonization spectrum. The general form of the 
RTE for the photon occupation number n{i') = J(i')/v^ 
with subrelativistic electron temperature in the presence 
of bulk motion was first derived by Blandford & Payne 
(1981, hereafter BP81). Later, Titarchuk et al. (1997, here- 
after TMK97) showed that analytical solutions can be 
found if the velocity profile (assuming spherical symme- 
try) of the matter follows the free-fall law, vr oc R~^/'^. 
In this case, the equation can be written as Lxn{x,T) + 
Lrn{x,T) = —s{x,t) and the solution can be obtained 
using the variable-separation method in the form 

oo 

n{x,T) ^^CkRk{T)Nk{x), (1) 
fc=i 

where Ck and i?fe(r) are the expansion coefBcients and 
eigenfunctions of the space operator Lt, respectively, 
while Nk{x) is the solution of the differential equation 

L:cNk{x) - ikNk{x) = -s{x), (2) 

where 7^ oc and is the k*^-eigenvalue of the space 
operator. 

The Comptonization spectrum is mostly dominated 
by the first eigenvalue (sec Fig. 3 in TMK97), while the 
terms N^ix) with k > 2 represent the fraction of pho- 
tons that leave the medium without appreciable energy- 
exchange. Starting from the results of TMK97, Farinelli 
et al. (2008, hereafter F08) developed a model (comptb) 
for the X-ray spectral fitting package xspec, which com- 
putes the emerging spectrum by means of a numerical 
convolution of the Green's function of the energy operator 
with a blackbody (BB)-like seed spectrum. The model has 
been successfully applied to a sample of low-mass X-ray 
binaries hosting a neutron star (NS). The method of the 
variable separation has also been adopted by Becker & 
Wolff (2007, hereafter BW07), to find analytical solutions 
of the RTE in the case of cylindrical accretion onto the 
polar cap of a magnetized NS. The starting equation in 
BW07 is formally the same as in BP81: the most signif- 
icant difference is that the Thomson cross-section is re- 
placed by an angle-averaged cross-section that takes into 
account the presence of the magnetic field (B ~ 10"'^^ 
G). Following the results of Lyubarskii & Sunyaev (1982), 
BW07 assumed a velocity profile v{t) oc — t, which al- 
lowed the RTE to be separable in energy and space. Note 
that the assumed velocity profile implies that the matter 
flow stagnates at the stellar surface, which is at odds with 
the solution of TMK97, where the matter velocity is in- 
creasing towards the central object, which can be either 
a NS or a black hole (BH). When the velocity profile is 



not free fall-like (TMK97, F08), or oc r (BW07), the vari- 
able separation method can no longer be applied and the 
solution of the RTE can be obtained only with numerical 
methods. 

We report a numerical algorithm that allows the solu- 
tion of the RTE in the FP regime using finite-differences 
for any desired velocity profile and seed photon spatial and 
energy distribution. We apply in particular it to cylin- 
drical accretion towards the polar caps of a magnetized 
NS, following the approach of BW07. The algorithm es- 
sentially uses a relaxation method, therefore it finds the 
asymptotic (stationary) solution of the RTE for a given 
initial value (at time t = 0) condition. Our work is struc- 
tured as follows: in Section 2 we describe the kernel of the 
algorithm for generic two- variable elliptic partial differen- 
tial equations; in Section 3 we formulate the problem for 
the general RTE and appropriate boundary conditions; in 
Section 4 we consider the more specific case of a system 
configuration with azimuthal symmetry, typical of cylin- 
drical accretion; in Section 5 we show the emerging spec- 
tra obtained for different sets of the theoretical parameter 
space; finally, in Section 6 we briefiy discuss possible as- 
trophysical consequences and implementations (e.g., for 
xspec) derived from the application of the algorithm. 



2. The general elliptic partial differential 
equations 

The algorithm we report is essentially based on the relax- 
ation methods, which allow one to find the solution of a 
boundary elliptical problem. The differential equation has 
to be written by finite differences. Once the sparse matrix 
is defined, it can be split into layers over which an itera- 
tion process is applied until a solution is found (Press et al. 
1992). The general form of a linear second-order elliptic 
equation with vanishing mixed derivatives and a source 
term can be written as 

^(^' y) ^ + y) ^ + y)u + W{x,y)^ 

„, .du du „, . 
+Z{x,y)— = — -S{x,y). (3) 

We define a three-dimensional grid of discrete points for 
the variables x, y and t\ 

Xi = xo+ihx, i = 0,1, . . . , Nx, 

yj = yo + jhy, j = 0,1,... ,Ny, 

tm = to + mht, m=l,2,...,M, (4) 

where h^.hy.ht arc the grid spacing. The function 
u(x, y, t) is evaluated at any point of the grid, so we write 
it as wf". We write the first and second derivatives over 
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the variables using finite differences: 



the system (9), which gives 



du 
dx 



dy 
du 

'at 



2u] 



hi 



u"" — u 



m— 1 



(10) 



lly 



j.m i.m — 1 



The third term on the right-hand side of equation (10) 
represents the residual error in the numerical solution. As 
(5) a first step, we collect the terms with the same index in 
both equations (9) and obtain 



Substituting the above definitions into equation (3) and 
collecting terms, we obtain 



i I l7 i-'in , 7 i,m , ij 7 — l,m 

cLiUi_-i + tr-ui + c^i^i+i + "i 

/ j.m 



1 7 7-^ 



j.m — l\ 

i ) 



q3 











ht) 



m-l/2 



V 



,,m — 1 



-1/2 



ht 



(11) 



where 

3 _ V{xi,y^) 

J _ V{xuyi) 
'~ hi ^ 

_ 2W{xi,y^) Z{x,,yi) 

hi ^ hy 

• _ >V(x„j/J) _ Z{x,,y^) 

"y "'2/ 



Both equations are defined inside a 2D (a;, y)-domain, 
(6)with boundary conditions defined according to the spe- 
cific problem under consideration. First, for any m and 
j values, we must impose the boundary condition on the 
left-hand side of the a;-domain {i = 0) for the function 



J,TO-l/2 



.J,m-l/2 



hx 

Qixi,y^) 
hx 



9o^ 



(12) 



while the source term Sq is defined at the beginning. Thus, 
for i = 0, equation (6) can be written as 



u: 



j,m-l/2 




j,m-l/2 



(13) 



where 



r-' — — 



ht 



-"■0 



J_ 

ht 



Si = S{xi,y^). 

The operators over 
as 

Ax = (4 + 14 + 4, 



(7) 



nj,m-l 



4 + ei + f^ + -)ur-'-Si, (14) 



J . , , J.1 J j2 J with the coefficients determined in equation (7). For 
the X and y variables are then defined . , . , , . 



using equation (13) we obtain 



di+ei + fi. 



(8) 



w 



i,m-l/2 



J, m-l/2 



The solution procedure consists of dividing equation (6) ^^-^^ ^^.^^^^ 



into two equations. The first gives us the solution for an 
iiiterniediate m — 1/2 layer, while the second provides the 
solution for the m-layer. As starting point we need to es- 
tablish an initial guess for the function at m = 1. The where 
system of equations to be solved is thus 

= 



^i,m-l/2^2iuj™-'/' + .ft'f, 

-^3 _ £l !-''■() 



Ki 



AxU 



m-l/2 



Ayiu"" - U""-^) = 



u^ - 



Iterating the process, we obtain the general form 



ht 



(9) 



7-,m-l/2 _ f J 7,m-l/2 r^j 



(15) 
(16) 

(17) 
(18) 



in which we have temporarily dropped the indices 
From the system we notice that the intermediate layer 
?Ti — 1/2 is needed only to build up the solution at the 
subsequent layer in m. The numerical accuracy of the so- 
lution can be estimated by combining both equations in 



where 

Li = 



„] TO 



1_ 

ht 



(19) 
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At the right boundary of the a;-domain {i = N^), we im- 
pose the second boundary condition 

4?"'^'= 5k- (20) 

Now, using equation (18) we can thus build up the solution 

over the x-variablc itcrativcly as 

j,m-l/2 



U 



j,m-l/2 
JVx-2 



-1/2 



r-l -I- K-f 



fj J".»n-l/2 
^N,-2"'N,-l 



^ JV.-2' 



(21) 



Therefore, the construction of the solution is obtained in 
two steps: a bottom-up process which allows one to build 
the coefficients Ll and Kl (Eq. [19]) starting from the left 
boundary condition on Uq (Eq. [12]), followed by a top- 
down procedure determined by the right boundary condi- 
tion (Eq. [20]). 

Once the solution over the a;- variable for the m — 1/2 
layer is obtained for any j (the index of the y variable), we 
then seek the solution of the second equation in the system 
(9) by following the same procedure described above with 
initial boundary condition at j = 

u^r = Lyr + Kl (22) 
and, similarly to equation (19) 



where 



fi 



ht 



di 



f + — 



j.m—1 



+ ey + 



^ j,m-l/2 



ht 



(23) 



(24) 



depends on the solutions u^''""^^^ and uf"^~^ obtained in 
the layers m — 1/2 and m — 1. 

As required for the procedure over the a;- variable, the coef- 
ficients and Kf , built from L° and Kf, arc determined 
by the left boundary condition (j = 0) for the function 
u?, and the solution for any j is determined by the right 
boundary condition uf^ = gf^ : 



-1,771 £j^^ 



Ny-2,m 



0,771 



-1 Ny,m . £>Ny- 

u- -I- K. " 



-2 Ny 



-hm^^Ny 



(25) 



After constructing the solution over the x and y variable, 
the solution of the top layer m becomes the initial function 
of the bottom layer related to iteration m + 1 according 
to the scheme 

TO = 1 ^ (?i°,ul/^ul), 

m = 2 (u\u3/^m2), 



It is also worth mentioning that at the first iteration m = 1 
an initial guess function u{'^ must be assumed, which is 
needed to find the solutions ic'i^^'^ and u{'^ in the system 
(9). The loop over to stops when the same convergence 
criterion is satisfied, which physically means that the so- 
lution has "relaxed" to its stationary value. One possible 
criterion could be 1 - £ < |u|'"~-^|/|u|™~-^/^ < 1 -|- £ and 
1 — £ < |m|™~-^/^|/|m|™ < 1 -h £, where £ is a user-defined 
numerical tolerance. In the next section we will show in- 
stead another convergence criterion we have chosen for 
stopping the iteration procedure, for the particular case 
of the RTE. 

3. Application to the radiative transfer equation 
and boundary conditions 

The general form of the RTE in the presence of subrela- 
tivistic bulk motion for a plasma with constant tempera- 
ture Te is given by (see Eq.[18] in BP81) 



dt 



Vn 



3ne<T(l') 



Vn 



1 — dn 



}_d_ 
dv 



me 



n - 



dn 



(27) 



where n(i/, r) is the zero-moment occupation number of 
the radiation field intensity, V is the plasma bulk veloc- 
ity vector, a{v) is the electron scattering cross-section, 
ne(r) is the electron density and j{v, r) is the source term. 
Because the spectral formation is determined by the op- 
tical depth r of the system, we use the latter quantity as 
the actual space variable. The solution of equation (27) is 
fulfilled by imposing the boundary condition at the sur- 
face defined by r = 0, (which represents the starting point 
of the integration domain) for the spectral flux, which is 
given by 



1 



3ncCr(i^) 



Vn 



1 dn 
— Vi'—— 

3 di^ 



(28) 



Under particular symmetries of the system configuration 
(e.g., cylindrical or spherical), the problem becomes one- 
dimensional. For constant electron temperature To it is 
also more convenient to use the adimensional variable 
X = hv/kTe, moreover, when performing numerical inte- 
gration using finitc-difFcrcncc methods, we use a logarith- 
mic binning of the energy through the additional change 
of variable a; ^ e*. Under these assumptions, equation 
(28) becomes 



Fiq,T) = 



\dJ_ 
3lh 



-V 



dJ_ 
dq 



J 



(29) 



where J = n is the specific intensity. 

At the inner boundary we impose the condition 



M-l „M-l/2 



(26) 



(30) 
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where A is the albedo at the surface. A fully absorptive 

surface (A = 0) is appropriate for a BH, while < A < 1 
accounts, e.g., for a NS atmosphere. However, the in- 
ner boundary condition (30) depends on energy as well 
as space (see Eqs. [28] and [29]). For mixed boundary 
value problems, no analytical solutions are possible (see 
Appendix E in TMK97) and numerical methods prove 
to be unstable. However, in the energy range where the 
spectrum is a powerlaw J{x,t) = R{t)x~°', equation (30) 
becomes 

-i5+,.(„ + 3)fi = -2(l^)fl. (31) 

where /3o is the bulk velocity at the inner radius (r = 0), 
and here the problem is reduced to a standard boundary 
condition over the space variable r. Writing the derivative 
in terms of finite-difference, equation (31) then becomes 

which can be written after collecting terms as 

= l + hr[l3o{a + 3) + G{A)]'''' ^^^^ 

where G{A) = 3/2(1 - A)/{l + A). We then set oui prob- 
lem as follows: first, because the solution of the RTE 
physically represents a specific intensity, it must by defi- 
nition be equal to zero in the limits E ^ and i? — > oo, 
therefore we set Uq = = (see Eqs. [12] and [20]). For 
the behavior of the function tij for r = (j = 0), equation 
(33) immediately allows us to define (see Eq. [22]) 

^ l + hr[l3o{a + 3) + G{A)y = °- ^^^^ 

As outer boundary condition over r, we impose that u^^ = 
0, which means that the specific intensity goes to zero for 

We emphasize that the condition > for any (i,j)- 
value implies a specific restriction in the choice of the step 
size hr, which ensures that L° > (as /3o < 0). More 
specifically, we imposed the condition on hr such that the 
number of steps over t he Nt = Tmax//ir > 10- 

3.1. The iteration procedure 

As already mentioned in Section 2, it is necessary to choose 
a convergence criterion for stopping the iteration over the 
m variable. We proceeded in the following way: at each 
iteration m, we computed the spectral index am of the 
solution u^'™ (corresponding to r = 0) in a given energy 
range E^i^i — -Emax- To minimize bias or wrong estimate 
of am, the definition of the energy interval for the com- 
putation of the spectral slope must be chosen carefully. If 
the seed photon spectrum is a BB with temperature fcTbb, 
a reasonable choice can be the assumption iSmin ~ 7A;Tbb 
and i^max ~ 20A;Tbb, respectively, given that this interval 



is above the major contribution of the BB component and 
below the expected high-energy cut-off value. 

Once am is estimated, it is inserted into equation (34), 
which accordingly represents the boundary condition at 
T = for the iteration m+1 (see Eq. [22]). We then com- 
puted the new index Um+i for u^''"'*'^, and again inserted 
it into equation (34) at iteration m + 2, and so on. The 
routine is stopped when a„i and am+i differ less than 
10~^ provided that the condition holds for a sufficiently 
high number of iterations (> 100). Note that the same 
criterion is adopted also if /3o = 0, even if then of course 
I/° remains constant across the iteration. We have also 
verified that this criterion automatically also satisfies the 
convergence of the norms juj™"^, \u\"^~^/^, and \u\"^. 

4. Cylindrical accretion onto a magnetized 
neutron star 

We applied our algorithm to solve the RTE for accretion 

towards the polar cap of a magnetized NS, whose mathe- 
matical formalism was developed by BW07 in the frame- 
work of the spectral formation of accretion-powered X-ray 
pulsars. The relatively strong magnetic field {B ^ lO-'^^G) 
of the NS is expected to channel the accretion fiow to- 
wards the polar caps, and for low values of the altitude 
above the star surface, the problem can be treated in a 
axis-symmetric approximation where the space variable is 
defined by the vertical coordinate Z. The magnetic field 
moreover forces the medium to become birefringent as 
the effect of vacuum polarization, and birefringence en- 
tails the formation of two linearly polarized modes (or- 
dinary and extraordinary) of the photons, each having a 
characteristic scattering cross-section. For ordinary mode 
photons with energy below the first cyclotron harmonic 
at Ec « 11.57 Bi2 keV (where = B/IO^^ q)^ BW07 
defined angle-averaged cross-sections parallel and perpen- 
dicular to the lines of the magnetic field as (T|| = 10~^crT 
and a± = (Tt, respectively, where (Tt is the Thomson scat- 
tering cross-section. This is indeed the only approximation 
that allows to treat the problem analytically or numeri- 
cally. We note that Ferrigno et al. (2009), starting from the 
analytical solutions reported in BW07, developed a model 
that was later almost successfully tested on the accret- 
ing pulsar 4U 0115+63. Their model is based essentially 
on the convolution of the column-integrated Green's func- 
tion of the thermal plus bulk scattering operator with a 
given seed photon distribution. The basic assumption of 
this derivation is that the velocity profile of the accreting 
matter is assumed to be v{t) oc — t, which allows one to 
find analytical solutions through the variable separation 
method (Eqs. [36] and [37] in BW07). The numerical al- 
gorithm we developed directly solves the RTE, without 
the need of this prescription for the dynamical configura- 
tion of the accreting matter field, and we included some 
modifications with respect to the approach of BW07 and 
Ferrigno et al. (2009). 

First, following TMK97, we include in equation (27) a 
second term in the thermal Comptonization operator that 
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accounts for the contribution of the bulk motion veloc- 
ity of electrons in addition to their thermal (MaxwcUian) 
component. With this prescription in mind, equation (27) 
becomes 



To solve equation (37), it is necessary to define the behav- 
ior of the velocity profile /3(r). We considered two possi- 
bilities: in the first one, we assumed a general form 



(42) 



Idn 
c'dt 

n 



S{e,Z)^- 
n^a 1 d 



V dr 



dv e dn 
dZSc'de 



£2 de 



cdZ 

n + {kn + my/3) 



_d_ 

dZ \3nf.a 
dn 

a7 



1 dn 



dZ 



where the normalization constant is defined = 
(}q{Zq/ZsY\ and /3o is the terminal velocity at the alti- 
tude Zq. 

(35) The continuity equation for the system here considered 
gives the electron number density 



where e = hv, a = 10 ^ctt, while tesc is the photon mean 
escape timescale (see Eq. [17] in BW07) 



n^a^rl 



(36) 



Now, using the relation dr = ne(T\\dZ and the logarithmic 
binning of the adimensional energy x = hf/kT^, equation 
(35) becomes 



1 dJ S{q,T) 



npaucH dt 



H 



1 + 



meviry 



3kT^ 



+ 



+ 



e« - 35 - 



SkT^ie" - 3 + S) - m^viry 



dJ_ 

dq 



J + 



3kT^ 

1 v{T)dJ 



3H dr'^ He dr ' 



where we have defined the quantities 

a kT^ 



and 



H 



5 = 



C7|| nfleC 

1 d^(r) 



(37) 



(38) 



(39) 



3H dr ' 

where /3(t) = v(t)/c, while the dimensionless parameter 
^ is given by (see Eq. [26] in BW07) 

15.8 ro 



m 



(40) 



Equation (37) is given in the general form (3) and for this 
particular case, we have 



V{r) = 1 
Q{t) = 



mou(T)2 

3fcre ' 

3A;re(e9 - 3 -h 5) - met;(r)2 



7^(r) = e« - 3(5 



3fcTe 



Z{r) = - 
S{q,T) = 



v{j)_ 
He ' 
S(.q,T) 
H ' 



(41) 



M 



n= = 



^p\p{Z)\cRV 



(43) 



where m = M/Me, is the mass accretion rate in Eddington 
units and Rq is the radius of the accretion column. 
We then define the adimensional quantities z and ro 
through the change of variables Z — )■ RgQmz and Rq — ?• 
RsQmro, where m = M/Mq, while Msq and i?s© are 
the Sun mass and Schwarzschild radius, respectively. The 
effective vertical optical depth of the accretion column is 
then given by 



iz)=J\a,dZ' = cJl-, 



r? + 1 



(44) 



where C = 2.2xl0~^, and zq is the vertical coordinate at 
the NS surface. 

Inverting relation (44), we also define the velocity profile 
of the accreting matter as a function of the optical depth 
r instead of the space variable z 



/3(t 



(45) 



As a second possibility, following BW07, we considered 
the velocity profile 

m = (46) 

where * = 0.67^/^0 (see Eq. [32] in BW07). 

Given that in our model the optical depth t represents 
one of the free parameters, once it is provided in input 
together with adimensional accretion column radius ro, 
the accretion rate rn must be first computed either from 
equation (44), if /3(r) is defined as in equation (45), or 
from equation (28) in BW07 if ^(r) belongs to equation 
(46). This step is necessary to determine the ^ parameter 
(Eq. [40]), and requires fixing the maximum altitude of the 
accretion column 2;max- We assumed Zmax = 22;o, and all 
emerging spectra (see next section) were computed with 
this choice. 

5. Results 

In this section we report some examples of the theoreti- 
cal spectra obtained by the numerical solution of equation 
(37) for different sets of the physical quantities that define 
the system. We consider a BB seed photon spectrum at 
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given temperature /cTbb with exponential spatial distribu- 
tion across the vertical direction, according to 

with the normalization constant defined as Cn = 
^km/^iO' where i?km and Dio arc the BB emitting area in 
kilometers and the source distance in units of 10 kpc, re- 
spectively. The spectra were computed using the velocity 
profiles defined in equations (45) and (46), respectively. 
The common parameters for both cases are consequently 
the BB temperature /cTbb; the electron temperature kTc, 
the optical depth r, the albedo at the inner surface A 
and the radius of the accreting column tq. On the other 
hand, for (3{t) belonging to equation (45), additional pa- 
rameters are the index 77 and the terminal velocity at the 
star surface /3o- We first present the results for this second 
physical case. 

In Fig. 1 we show the emerging spectra for different 
values of the electron temperature kTg and two terminal 
velocities (3o = 0.1 and /3o = 0.64. As expected, both times 
higher values of kTc produce flatter spectra and push the 
cut-off energy Ec to higher energies; on the other hand, the 
bulk contribution as a second channel of Comptonization 
depends on the value of kTc ■ The two extreme temperature 
values reported here, kTc = 5 kcV and kTc = 50 keV, 
are particularly instructive: for low electron temperatures 
the spectrum changes from BB-like when l3o = 0.1 to a 
cut-off power law with Ec 30 keV when /3o — 0.64, 
while the spectral change is much less enhanced for a hot 
plasma. These can be considered as typical examples of 
bulk-dominated and thermal-dominated Comptonization 
spectra, respectively. 

Together with the electron temperature, the optical 
depth T is an important parameter that plays a key role 
in determining the spectral slope and cut-off energy, as 
clearly shown in Fig. 2. We note that in Fig. 1 and Fig. 2 
the index of the velocity profile was chosen to be ry = 0.5, 
typical of accretion onto a compact object where gravity 
and radiation pressure are the only force terms that de- 
termine the dynamical configuration. Here, the terminal 
value of the matter velocity /3o depends on the ratio of 
the radiative and gravitational forces, provided the con- 
dition |Fr|/|Fg| <i 1 is satisfied. This relatively simple 
approach is valid for low values of optical depth r, while 
when r > 1 radiative transfer becomes important and the 
problem requires in principle a more accurate radiation- 
hydrodynamics treatment. 

It is outside the scope of this paper to compute the 
exact velocity profile for accreting matter under the pres- 
ence of a strong radiation field in a high optical depth 
environment. We merely introduced a simple parametriza- 
tion for modifying the velocity field by changing the index 
r], with the results shown in Fig. 3, for two different val- 
ues of the electron temperature kTc- As Fig. 3 shows, the 
lower the value of 77, the harder the spectrum: this behav- 
ior can be explained in a quantitative and a qualitative 



way. Indeed, as ry increases, the velocity profile (5{z) be- 
comes sharper, and for a fixed terminal velocity /3o, elec- 
tron temperature A;Te, and optical depth r, while photons 
diffuse through the bounded medium, on average the en- 
ergy of the electrons (caused by their Maxwellian plus 
bulk motion) will be lower, and consequently the net en- 
ergy gain of the photons due to inverse Compton will be 
less. From the mathematical point of view, it is worth 
mentioning that Mastichiadis & Kylafis (1992, hereafter 
MK92) reported the analytical solution of the RTE in 
the Fokkcr-Planck approximation with the variable sep- 
aration method for spherical accretion without magnetic 
field in the limit Te = 0. Assuming a general velocity pro- 
file oc , the authors showed that the spectral index 
of the fc"^-Comptonization order emerging spectrum yields 
Oik = 3 3Ak/(2 - ri) (see Eq. [1]), where Ak is the fc*^- 
eigenvaluc of the space operator. Using equation (B12) of 
MK92, it follows immediately that as 77 increases, the spec- 
tral index ak increases as well. This mathematical result 
in terms of spectral formation can be considered as general 
in the framework of the FP treatment, and is accordingly 
qualitatively meaningful for our research. 

Wc also emphasize that analytical solutions for rj ^ 0.5 
have been possible for MK92 only because of the condition 
Te = 0, which drops the thermal Comptonization operator 
in the RTE, while when Te > this is possible only for 
77 = 0.5 (TMK97, F08). 

In Fig. 4 we show results for different terminal bulk ve- 
locities ,5o for two electron temperature values. The figure 
can be considered an extension and completion of Fig. 1 
because more values of ^0 are shown to better appreciate 
the induced changes in the emerging spectra. 

The spectral modifications as a result of different val- 
ues of the albedo A at the inner surface are instead shown 
in Fig. 5, where we explored full absorption (^4 = 0) and 
full refiection {A = 1), together with other intermediate 
values. In the framework of a physical link to astrophysical 
objects it would be natural to associate a BH to the condi- 
tion A = and a NS to the condition A = 1, respectively 
(Titarchuk & Fiorito 2004; Farinelli & Titarchuk 2011), 
even though this latter assumption may be considered 
an oversimplification of the problem. A most realistic ap- 
proach would consist indeed in an energy-dependent treat- 
ment of the albedo, a problem that could be faced only 
with Montecarlo simulations, with the additional compli- 
cations arising from a detailed treatment of the star pho- 
tosphere (surface) properties. 

For our unavoidably simplified assumptions, the net 
effect of increasing values of A is a progressive flattening of 
the emerging spectra. This is physically explained because 
when A > 0, a fraction of photons (which becomes 100% 
when A=l) suffers on average more scattering with respect 
to A = 0. Qualitatively, the spectral modification leads in 
the same direction as an enhanced optical depth of the 
system. 

The last parameter that strongly influences the spec- 
tral formation is the radius of the accretion column ro, 
whose effects are shown in Fig. 6. Indeed, following the 
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BW07 prescription, the mean escape time for photons 
using the diffusion approximation t^sc oc Tq (sec Eq. 
[36]). On the other hand, both the bulk and thermal 
Comptonization parameters (j/buik and j/th, respectively) 
are related to the mean number of scatterings that pho- 
tons experience in the medium via 

?ybuik « A^^v""'Cbuik, (48) 
2/th ~ K%h, 

where A?^^""^, Cbuik, -^Iv and Cth are the averaged number 
of scatterings and the fraction energy gain per scattering 
for bulk and thermal Comptonization, respectively. Both 
-/V^""^ and are of course also proportional to tesc (see 
Eqs.[94]-[97] in BW07). 

Evidently therefore, for fixed velocity profile parame- 
ters £/ and T] (see Eq. [42]), once the optical depth r is de- 
fined (see Eq. [44]), to keep its value constant for increas- 
ing ro (as reported in Fig. 6), the accretion rate m must 
also increase in a way to keep the ratio m/rQ constant. 
Combining equation (36) and (43) yields tesc oc rh, which 
in turn leads to an enhancement of the Comptonization 
parameters ybuik and yth in equation (49) with a harden- 
ing of the spectral shape. 

Considering now the velocity profile defined in equa- 
tion (46) , we see that the results are qualitatively the same 
as in equation (45) as far as the spectral modifications in- 
duced by variations of kT^ are concerned (Fig. 7), r (Fig. 
8) and A (Fig. 9), respectively. But there arc opposite ef- 
fects that are induced in the emerging spectra by different 
values of the accretion column radius ro for the velocity 
profile here considered. Indeed, using equation (40) and 
the definition of r in equation (28) of BW07, which al- 
lows us to express the accreting matter velocity in terms 
of the z-coordinate, in spite of the optical depth r, it 

— 1/2 

is straightforward to see that /3{z) oc rg . In particu- 
lar, if zo = 2.42 and Zmax = ^zq we have /3max=0.60 for 
ro = 0.1, /3max=0.38 for m = 0.25, /3„,ax=0.27 for m = 0.5 
and /3max=0.2 for ro = 1, respectively. Note that because 
equation (46) describes matter that stagnates at the star 
surface, here 0rnax represents the velocity at the accretion 
column altitude Zmax- In other words, while using equa- 
tion (45), the choice of ro does not modify the velocity 
field of the accreting matter, which is only determined by 
the choice of /3o and ry, for (46) as ro increases the bulk 
contribution to the spectral formation becomes less impor- 
tant, and this drop is not compensated for by the increase 
of the photon mean escape time teso which, as explained 
above, would instead contribute to spectral hardening. 

6. Implementation in the xspec package 

Our model will be publically available and distributed as 
a contributed model to the official xspec^ web page. 
In Table 1 we report a summary of the free parameters of 
the model with their physical meaning. The code is written 

^ http://heasarc.nasa.gov/xanadu/xspec/newmodels.html 



Table 1. Parameter description of the xspec model compmag. 



Parameter Units Description 

T 

V 

/3o 

ro 
A 

Flag 
Norm 



using C-language, and can be easily installed following the 
standard procedure reported in the official XSPEC manual 
and in the brief cookbook, which will be delivered together 
with the source code. As a general concern for users, we 
point out that usually the emerging spectrum obtained 
from the Comptonization of a seed photon population with 
any given energy distribution S(E) can be presented as 
the sum of the seed spectrum and its convolution with 
the scattering Green's function G{E,Eq) of the electron 
plasma, each with their relative weight, according to the 
general formalism 

^(^) = -Xfl^^^^^^'^'' SiE)*G{E,Eo)], (49) 

where C„ is a normalization constant. The ratio A/ (A + l) 
is the Comptonization fraction, and its value qualitatively 
determines the contribution to the total spectrum of the 
Comptonized photons. The value of A, here not to be con- 
fused with the albedo, may depend on several geomet- 
rical and physical factors, such as the spatial seed pho- 
ton distribution inside the system configuration (see Fig. 
4 in TMK97). The lower the value of A, the more en- 
hanced the direct seed photon spectrum S(E). Examples 
of XSPEC models that use the definition in equation (49) 
are BMC (TMK97) and COMPTB (F08). Either model, how- 
ever, does not solve the full RTE including the photon 
spatial diffusion and distribution, the latter of which is an 
unknown quantity that is phenomenologically described 
through the continuum parameter log(A). On the other 
hand, it is not possible to change the value of log(A) ar- 
bitrarily in our present model, i.e., according to the ob- 
served spectra. Its value is implicitly determined once the 
seed photon spatial distribution is fixed. 

We presented the results of simulated spectra assum- 
ing an exponential distribution over r for S(E), which 
was assumed to be a BB; then, the transition from the 
low-energy part of the spectrum (the Rayleigh regime for 
E ^ SfcTbb) to the high-energy (Comptonized) powerlaw 
shape is almost smooth, which corresponds approximately 
to A ;:g> 1 in equation (49). Other seed photons spatial 



(keV) Seed photon blackbody temperature 
(keV) Electron temperature 

Optical depth of the accretion column 
Index of the velocity profile (Eq. [45]) 
Terminal velocity at the NS surface 
(Eq. [45]) 

Radius of the accretion column in 
units of the NS Schwarzschild radius 
Albedo at the NS surface 
= 1, /3(t) from equation (45) 
= 2, /3(t) from equation (46) 

-Rkm/-ClO 
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distributions can produce a diflFerent onset between the 
BB peak and the powcrlaw-likc regime. In general, for 
observed spectra where a direct and enhanced BB-like 
component is required by the fit, our claim is to model 
the source continuum with modclization of the type BB + 
COMPMAG by preferably keeping equal to each other the 
temperatures of the direct and Comptonized BB compo- 
nent. 

7. Conclusions 

We have developed a numerical code for solving elliptic 
partial differential equations based on a relaxation method 
with finite differences. In particular, we reported a spe- 
cific application of the algorithm to the radiative transfer 
equation in the Fokker-Planck approximation, which is of 
particular interest for high-energy astrophysical applica- 
tions. We considered cylindrical accretion onto the polar 
cap of a magnetized neutron star, using the mathemati- 
cal formulation of the problem reported in Becker & Wolff 
(2007) with some modifications. Specifically, we included 
the second-order bulk Comptonization term in the scatter- 
ing operator and we considered different velocity profiles 
for the accreting matter. 

The code for the computation of the emerging spec- 
tra in this configuration has been written with the aim 
to implement it in the xspec package and will be deliv- 
ered to the scientific community. Because angle-averaged 
cross-sections caused by the magnetic field were included, 
the model is suitable to be applied to the observed spectra 
of sources where most of spectral formation is claimed to 
form close to a NS with high magnetic field {B ^ IQ-'^ G), 
such as accreting X-ray pulsars and supergiant fast X-ray 
transients. Of course there are some unavoidable simplifi- 
cations in the model, such as the assumption of constant 
electron temperature of the accretion column. We note, 
however, that the isothermal condition is typical of any 
Comptonization model implemented in XSPEC because it 
is fundamental for users to reach a compromise between 
the accuracy of the physical treatment and the computa- 
tional speed. More accurate theoretical investigations of 
the accretion processes can be performed only with MHD 
simulations performed through PC cluster-computing re- 
sources, which are beyond the scope of the present work. 

As in many theoretical models, the number of available 
free parameters is higher than the number of the observ- 
able ones. Therefore a correct working approach is to keep 
some parameters fixed during the X-ray spectral fitting 
procedure to avoid degeneracy. 

If this model is used in any publication, we kindly ask 
the authors to cite this paper in the reference list. 
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Energy (keV) Energy (keV) 

Fig. 1. Emerging spectra obtained from the solution of equation (37) for different values of the electron temperature kTc and 
velocity profile defined in equation (45). In both cases the fixed parameters are fcTbb = 1 keV, r = 0.2, rj — 0.5, ro = 0.25, 
A — 1. Left panel: /3o = 0.1. Right panel: /?o = 0.64. 




Fig. 2. Same as Fig. 1 but for different values of the optical depth r. Fixed parameters are kThh ~ 1 keV, rj = 0.5, /3o ~ 0.64, 
ro = 0.25, A = l. Left panel: kTc = 5 keV. Right panel: kTc = 15 keV. 




Fig. 3. Same as Fig. 1 but for different values of the index of the velocity profile. Fixed parameters are kThh = 1 keV, r — 0.2, 
/3o = 0.64, ro = 0.25, A^l. Left panel: kTc = 5 keV. Right panel: kTc = 15 keV. 
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Fig. 4. Same as Fig. 1 but for different values of the inner velocity Po- Fixed parameters are kThb = 1 keV, r = 0.2, 77 — 0.5, 
ro = 0.25, A^l. Left panel: kT^ = 5 keV. Right panel: kT^ = 15 keV. 
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Fig. 5. Same as Fig. 1 but for different values of the albedo A, with the velocity profile of equation (45). Fixed parameters are 
fcTbb = 1 keV, r = 0.4, r/ = 0.5, /3o = 0.64, ro = 0.25. Left panel: kT^ = 5 keV, Right panel: kT^ = 15 keV. 
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Fig. 6. Same as Fig. 1 but for different values of the accretion column radius ro, with the velocity profile of equation (45). Fixed 
parameters are kThh = 1 keV, r = 0.2, r) — 0.5, /3o ~ 0.64, A — 1. Left panel: kTc — 5 keV. Right panel: kTe = 15 keV. 
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Fig. 7. Emerging spectra obtained from the solution of equation (37) for different values of the electron temperature fcTc, with 
the velocity profile of equation (46). In both cases the fixed parameters are fcTbb = 1 keV, ro = 0.25, A = 1. Left panel: r = 0.2. 
Right panel: r — 0.4. 
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Fig. 8. Same as Fig. 1 but for different values of the optical depth r, with the velocity profile of equation (46). Fixed parameters 
are fcTbb = 1 keV, ro = 0.25, and A= 1. Left panel: kTc — 5 keV. Right panel: fcTc — 15 keV. 
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Fig. 9. Same as Fig. 1 but for different values of the albedo A, with the velocity profile of equation (46). Fixed parameters are 
fcTbb = 1 keV, r = 0.4, and ro = 0.25. Left panel: kT^ = 5 keV. Right panel: kT^ = 15 keV. 
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